Mean Field Theory of Polynuclear Surface Growth 
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We study statistical properties of a continuum model of polynuclear surface growth on an infinite 
substrate. We develop a self-consistent mean-field theory which is solved to deduce the growth 
velocity and the extremal behavior of the coverage. Numerical simulations show that this theory 
gives an improved approximation for the coverage compare to the previous linear recursion relations 
approach. Furthermore, these two approximations provide useful upper and lower bounds for a 
number of characteristics including the coverage, growth velocity, and the roughness exponent. 
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I. INTRODUCTION 

Kinetics of surface growth is a fascinating field that 
has been the subject of intense current research 
It is well established that as the surface grows its mor- 
phology remains scale invariant, and for example, fluctu- 
ations in the interface height exhibit an asymptotic scal- 
ing behavior. While the understanding of growth on one- 
dimensional (Id) substrates is rather comprehensive, in 
the physical case of two-dimensional substrates current 
theoretical understanding remains incomplete . 

In this study we focus on an appealingly simple yet 
non-trivial surface growth problem, the so-called polynu- 
clear growth (PNG) model §-§. The PNG model de- 
scribes the evolution of islands that nucleate at random 
on top of previously nucleated islands and grow in the ra- 
dial direction. This model is appropriate for describing 
situations where there is a competition between growth 
along the step edges and growth due to nucleation, as is 
the case in polymer crystal growth j^jjj . 

The submonolayer version of the PNG model is iden- 
tical to the classical Kolmogorov-Avrami-Johnson-Mehl 
(KAJM) nucleation-and-growth process flf]|-|l^l where 
exact results for the coverage and more detailed statisti- 
cal properties are possible p3| |l6| . Using the exact solu- 
tion for the KAJM coverage, an approximate Linear Re- 
cursion Relations (LRR) approach to the PNG process in 
arbitrary dimension was suggested P,p7|-^9[. Addition- 
ally, steady states on Id substrates were obtained ana- 
lytically using the fact that kinks (antikinks) are uncor- 
rected p, p0pi[ |. However, the non-equilibrium behavior 
and especially the asymptotic time dependence remains 
an open problem, despite a number of studies [i"g|-|24|]. 

Our goal is to develop a self-consistent approach that 
can be viewed as a mean field theory (MFT) of the PNG 
process, and to compare it with the previous LRR ap- 
proximation as well as with numerical simulations. We 
will show that MFT offers a better description for the 
time dependent coverage, and that the two approaches, 
when combined, provide upper and lower bounds for the 
coverage and the growth velocity. 

The rest of this paper is organized as follows. In the 



next section we define the PNG model. In Sec. Ill, we de- 
velop a self-consistent mean field approximation. We find 
that possible growth velocities are bounded from below, 
v > fmin, and show that the minimal velocity is actually 
selected. We also solve for the coverage profile in the 
tail regions. In Sec. IV, we review the LRR approach 
and derive the coverage profile analytically. In Sec. V, 
we present simulation results. Conclusions are given in 
Sec. VI. 



II. THE PNG MODEL 

The PNG model is defined as follows. Consider a flat 
uniform G?-dimensional substrate at time t — 0. Seeds of 
negligible size nucleate randomly at a constant rate per 
unit area, and grow with a constant velocity in the ra- 
dial direction. When two islands on the same layer meet 
they coalesce, and the joint perimeter continues growing 
in the corresponding radial direction. Meanwhile, nucle- 
ation continuously generates additional layers on top of 
previously nucleated layers. Clearly, there are no over- 
hangs in this model, a feature that considerably simpli- 
fies the analysis. Another important simplification in the 
PNG model is that the nucleation rate is uniform in time 
as well as in space, i.e., it is independent of the local 
surface structure. Without loss of generality, we set the 
nucleation rate and the radial growth velocity to unity. 
This can be achieved by an appropriate rescaling of space 
and time. 

In this study, we concentrate on Sj(t), the uncovered 
fraction in the jth layer at time t. This important char- 
acteristic of multilayer growth gives the net exposed frac- 
tion of the jth layer, Sj(t) — Sj-\(t), and therefore can 
be used to calculate relevant statistical properties. In 
general (/(j)) = £~ f{j)[S j+1 {t) - Sj(t)], and in par- 
ticular the average height is given by 

oo 

h(t) = {j) = J2j[S j+1 (t)-S j (t)]- W 
j=i 

Fluctuations in the height are quantified by the mean 
square width or roughness, w 2 (t), 
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w 2 (t) = if) - or 



(2) 



We expect a linear growth in time for the average height, 
h(t) ~ vt, and an algebraic growth for the interface 
width, w(t) t" , with a priori unknown roughness ex- 
ponent p. In other words, the uncovered fraction obeys 
the following wave-like form 



S j (t) = F 



j - vt 



(3) 



The argument reflects the overall shift in the position of 
the wave-front, and the multiplicative scale accounts for 
the algebraic widening of the front. 

Far from the front region, we anticipate the following 
extremal behavior of the scaling function F(z): 



1 — exp(- 
exp(— \z\ c 



oo; 
— oo . 



(4) 



The exponents a± thus characterize relaxation away from 
the front region. The exponent <r+ which describes large 
positive fluctuations in the height can be simply re- 
lated to the roughness exponent. Consider a large pos- 
itive height fluctuation, j = Avt, with A > 1. Such 
large "towers" can be created only by an anomalously 
large number, Avt, of nucleation events localized in the 
same region. Given the Poisson nature of the nucleation 
events, such fluctuations are suppressed exponentially. 
Thus, the quantity 1 — SAvt{t) is estimated by exp(— t), 
but since Eq. (Q) gives exp [— i <T +( 1- ' 3 )] , we conclude that 
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(5) 



III. MEAN FIELD THEORY 

In the following, we explore the non- equilibrium 
regime, i.e., we consider polynuclear growth on an in- 
finite substrate. Non-equilibrium behavior should agree 
with the early time behavior on finite substrates, since 
then finite size effects are still negligible. 



A. One Dimension 

We start with the PNG model in one dimension where 
a more comprehensive analysis is possible. In this situ- 
ation, steps nucleate in pairs and move away from each 
other with a constant velocity. The constant nucleation 
rate and the growth velocity are set to unity, without loss 
of generality. Hence, the length of an island at time t af- 
ter birth equals 2t. Consider fj(x,t), the density of gaps 



of length x at time t in the jth layer. This distribution 
evolves according to 



dfj(x,t) _ 2 dfj(x,t) 



dt 



dx 



(6) 



7i(*) 



-xf j (x,t)+2 dyfrfat) 



The spatial derivative term describes shrinkage of gaps. 
The last two terms account for changes due to nucleation 
and thus, are proportional to the overall nucleation rate 
at the jth layer, -fj(t). The loss term is proportional to 
the gap length and the gain term describes creation of 
gaps from larger gaps. 

Eqs. (|^) contain yet unknown nucleation rates Jj(t) 
which will be chosen to satisfy the correct kinetic equa- 
tions for the uncovered fractions 



dx x fj(x,t), 



and the gap (or island) densities 



Nj(t) = / dxfj(x,t). 



(7) 



(8) 



The uncovered fraction decreases with a rate propor- 
tional to the island density, Sj(t) — —2Nj(t), and by 
integrating Eqs. (^) we indeed recover this exact equa- 
tion. The island density changes due to disappearance 
of gaps as well as due to nucleation. The total nucle- 
ation rate is proportional to the exposed fraction of the 
jth layer and thus, Nj(t) = -2/j(0, t) + Sj(t) - Sj-i(t) 
On the other hand, 
AT i (t) = -2/ i (0,t) + 



by integrating Eqs. (|6|) we obtain 
jj(t)Sj(t). Therefore, the choice 



Sj-i(t) 



(9) 



guarantees that the gap density evolves according to ex- 
act rate equation. The expression of Eq. (^|) for the nu- 
cleation rate (t) in the jth layer is intuitively appealing 
since the total unit nucleation rate should be reduced to 
account for nucleation events below the jth layer. Taking 
into account that So (t) = we find 71 (t) = 1 and we no- 
tice that Eq. (||) for the first layer agrees with the exact 
KAJM equation |pj|| . Thus, the set of rate equations (||) 
with the nucleation rates (^) provides a self-consistent 
description of the Id PNG model. It is exact for the first 
layer, and additionally, the first two moments of the gap 
density satisfy the correct rate equations. However, it 
is a mean-field description since it assumes a spatially 
homogeneous nucleation rate jj (t) . 

The gap density is found to be exponential, and the 
formal solution reads 



fj(x,t) = g 2 At) exp 



-g j {t)x-2 / dTg 3 ( T ) 
Jo 



> (10) 



with gj (t) = f Q dr jj (r) . The uncovered fraction and the 
island density are evaluated using Eqs. (0) and (j|): 
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S j (t) = e- 2 fo dT3 ^ T \ 
Nj(t) = 9j (t)S 3 (t). 



(11) 



Evaluating d 2 In Sj(t)/dt 2 together with Eq. (||) and 
(jj{t) = Jj(t), leads to an infinite set of recursive dif- 
ferential equations for the uncovered fraction 



5', 



2(S i -5 J _i) = 0- 



(12) 



Eqs. ( fL2| ) should be solved subject to the initial condi- 
tions Sj(Q) — 1 and Sj(Q) = for j > 1. The recursive 
structure of Eqs. ( |l2| ) reflects the fact that kinetics of a 
given layer is unrelated to that of all layers above (thanks 
to the absence of overhangs). Additionally, Eqs. jl^ ) 
agree with the nature of the PNG model implying that 
kinetics of a given layer should be directly coupled only 
to the previous layer. 

Using So(t) = 0, Si(t) is determined, then 5*2 (i), etc. 
Of course, for the first layer, the KAJM nucleation-and- 
growth results are reproduced |13| 



Si(t) = e~ r , 

h{x,t) = ee- xt - t2 



(13) 



It is also possible to solve analytically for the second layer 



S 2 (t) = cosh 2 ie~* 



(14) 



f 2 (x, t) = [t cosh t - sinh t} 2 e -^- tanh ^ x ' t2 



Using the transformation, Sj(t) = exp \uj(t) — i 2 ] , the 
differential equations (12) formally simplify to a directed 
version oftheToda equations p5| , ilj = 2exp[v,j—i — Uj], 
with the initial conditions Uj{0) = %(0) = and the 
boundary condition u\{t) — 0. Despite this simplifica- 
tion it is not possible to integrate these equations, and 
we solve numerically for Sjif). Fig. 1 shows how the cov- 
erage in a given layer changes with time. We see that the 
coverage quickly relaxes onto a traveling wave form with 
a finite width, Sj(t) — > F(j — vt). 

Some quantitative characteristics of the traveling wave 
can be determined analytically. For j — vt 3> 1, the non- 
linear term in ( p"2j ) is negligible and Eqs. (]l^) become 
linear. Thus, we write 



l-S,-(t)~e 



-<x{j-vt) 



j-Vt> 1, 



(15) 



with a yet unknown coefficient a. 
Eqs. (Ill) gives 



Substituting into 



(16) 



The right hand side has a minimum at a = 1.59362. 
Therefore any velocity in the interval [i; m ; n ,oo) with 
^min = 1-75735 is possible. Our numerical integration 
shows a velocity that falls within 0.1% of v m in, thereby 
implying that the minimal velocity is indeed selected. 



Such minimum velocity selection is ubiquitous and oc- 
curs for a wide class of initial conditions [EfAPTf . 
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Fig.l The uncovered fraction Sj(t) vs. time for layers j — 
20, 40, 60, and 80. Clearly, the coverage follows a traveling 
wave form. 

To obtain the asymptotic behavior in the other ex- 
treme, z = j—vt — > — oo, we first note that ^(t) ^ Si(t), 
as follows from ( |T3| ) and (14). We assume this for all lay- 
ers far behind the front; we will check this assumption 
Sj(t) ^> Sj-i(t) a posteriori. This reduces Eqs. ( p^ ) 
to Sj - Sl-Sj 1 + 2Sj = which is solved to yield 
Sj(t) — exp(— t 2 + Ajt + Bj). The traveling wave form 
implies that Sj(t) should be a function of a single vari- 



able z 



vt. This determines the constants Aj and 



Bj, and we find Sj(t) — F(z) ~ exp (— z 2 /v 2 ). One can 
verify that the assumption Sj(t) 3> Sj-\(t) is valid. 



B. Arbitrary Dimension 

The above analysis cannot be generalized in a straight- 
forward manner to d ^ 1 since the gap distribution is 
intrinsically one-dimensional. However, it is still possi- 
ble to obtain a mean field description for the uncovered 
fraction. 

Consider first a KAJM nucleation-and-growth process 
where the rate of nucleation events 7(i) is homogeneous 
in space but time dependent. Ignoring overlap between 
growing disks, the uncovered fraction S(t) decreases with 
time according to 



dS(t) 
dt 



= -dfl d / dr 7 (r)(< -r) 



d-l 



(17) 



where fl d = 7r d / 2 /r(l + d/2) is the volume of the d- 
dimensional unit sphere and dfl^ is its surface area |2§| ] 
(r is the gamma function). Of course, Eq. ( p7| ) overes- 
timates the decay rate since some of the area is already 
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covered. Nevertheless, this may be corrected if (—dS/dt) 
is reduced by the uncovered fraction, S, 



dfl d S(t) / dr 7 (r)(<-r) 



dS{t) 
dt 

For a constant nucleation rate, 7 = 1, Eq. 

n d t d + v 



d-l 



(18) 



gives 



S(t) = cxp 



1 



(19) 



Thus, the exact KAJM coverage is recovered. Further- 
more, a generalization of the KAJM solution to time- 
dependent nucleation rates appears to be equivalent to 
Eq. @ Q. 

We now return to the PNG model. Note first that 
the nonlocal in time integro-differential Eq. (|l8|) can be 
converted to a higher order ordinary differential equa- 
tion d d+1 \nS 3 {t)/dt d+l = -d\n d ^{t). Unlike the Id 
case, it is not possible to derive the nucleation rate self- 
consistently. However, assuming a spatially homoge- 
neous nucleation rate |tJ implies the total nucleation rate 
Jj{t)Sj{t) = Sj(t) - £j_i(t) and therefore Eq. (|). We 
thus arrive at the following generalization of the mean- 
field Eq. (F2I) for the uncovered fraction 



d d+\ 

dt d+1 3 



d\ CL, 



1 



Sj-i(t) 



0. 



(20) 



The analysis presented in the one-dimensional case ap- 
plies for arbitrary dimensions. For example, the trans- 
formation Sj(t) = exp [uj(t) — Qdt d+1 /(d + f)] reduces 
Eqs. ( ^o"l) to a set of generalized directed Toda equations 
d d+1 Uj /dt d+1 — dlQd exp[uj^i — Uj]. Analysis of these 
equations or Eqs. (|20| ) reveals that the coverage relaxes 
to a traveling wave with a finite width. To determine 
the growth velocity, we insert the ansatz of Eq. (|l5|) into 
Eqs. <M) to find 



,,d+l 



d\Q d 



d+l 



(21) 



This provides the lower bound for the growth veloc- 
ity Vd > v d lln - Again the minimal velocity should 
be selected, and for exam ple vq = 1, v% = 1.75735, 
V2 = 1.67115, and Vd — \j2ne/d when d — > 00. Fur- 
thermore, in finite dimensions mean field theory predicts 
universal exponents (3 — 0, au. = 1, and er_ = 2. (Note 
that the relationship of Eq. (||) is obeyed) . 

In the zero-dimension limit, the behavior changes qual- 
itatively. Indeed, for the d — case Eqs. (^o|) become 
linear, 



dSj_ 
dt 



S 7" — >S* <i 



(22) 



Solving ( P4 ) recursively yields the uncovered fraction, 
Sj(t) = e _ *J^T £*/«!. Alternatively, by treating the 
variable j as continuous, this set of linear equations re- 
duces to a simple convection-diffusion equation with a 
unit velocity and a diffusion coefficient D = 1/2. Conse- 
quently, the roughness becomes diffusive, i.e., (3 = 1/2. 



IV. LINEAR RECURSION RELATIONS 



The Linear Recursion Relations (LRR) approach em- 
ploys the exact uncovered fraction S\(t) in the first layer, 
provided by the KAJM solution © Nu- 
cleation in the (j + l)st layer proceeds only on the al- 
ready covered fraction of the jth layer, formed with 
rate —dSj/dt. Subsequent covering proceeds as in the 
KAJM, so one anticipates that nucleation events in the 
time interval (r, r + dr) make a contribution S\ (t — 



-Sj(r) dr to the exposed fraction Sj + i(t) — Sj(t) 

in the (j + l)st layer. This leads to a recursion relation 
between adjacent layers: 



S j+1 (t) =Sj{t)- f drS^t-r) 
Jo 



dr ' 



(23) 



The first layer coverage ( |l9| ) can be recovered by setting 
the substrate coverage appropriately, dSo/dt = —d(t). 

Although MFT and LRR are both recursive as every 
layer is coupled to the preceding layer, they differ in that 
the MFT equations are nonlinear, while the LRR equa- 
tions are linear. Nevertheless, when d — > 0, both approx- 
imations are identical. Indeed, multiplying Eq. (|2^) by e* 
and differentiating, one recovers the MFT equation (|2^) . 

Using the Laplace transform, analytical results for the 
growth velocity and the interface width have been estab- 
lished |3(^,|l8],|l9| . Below, we give an alternative and sim- 
pler derivation which additionally provides the asymp- 
totic behavior of the coverage profile. In the long time 
limit, it is reasonable to treat the layer number j as a con- 
tinuous variable. Replacing the difference Sj+i — Sj by 
a partial derivative, the recursive relations (j23|) become 



dS_ 



'drS^r) 98 ^ 



dt' 



(24) 



1_55 

v d dt 



j_ 1 y^+fj d 2 S 

2u| T 2 (d+2\ ~dj?' 



with Vd the growth velocity 



Vd 



a, 



\d+l J 



d + 2 



(25) 



The second line in ( p4| ) has been derived by expanding 
S(j, t — t) in a Taylor series in r, keeping only the two 
dominant terms of the expansion, and replacing the up- 
per limit in the integral by 00. The following growth 
velocities are found: vq = 1, v\ = 2/y / 7r = 1.12838, 
v 2 = 1.13719, and v d — > \j2-nejd as d ~ * 00. The veloc- 
ity is almost constant for physical dimensions (it varies 
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by less than 4% in the range 1 < d < 4) indicating the 
weak dimension dependence of this approach. 

Changing variables from ( j, t) to ( j, £ = j — Vdt) recasts 
Eq. into a diffusion equation 



dS_ 
~dt 



D 



8 2 S 



(26) 



The constant D = (|±f) jv 2 (j±f) plays the role 

of a diffusion coefficient and controls the width of the 
interface. In obtaining Eq. (|26]), j was replaced by Vdt. 
This is clearly valid in the scaling limit, j — > oo, |£| — ► oo, 
j ~ £ 2 . The initial profile of the uncovered fraction, 
S(j,0), is a step function: SOLO) — for j < and 
S^j, 0) = 1 for j > 0. Solving (J2q) subject to these initial 
conditions yields 



S(j,i) = iErfc(- 



(27) 



with Erfc(z) — due u the error function p8[ . 

While the two approximations generally differ in their 
asymptotic behavior, they do agree in the extreme cases 
of d = and d = oo. 

In summary, the LRR approach predicts a dimension- 
independent "diffusive" width exponent (3 — 1/2. The 
shape of the coverage profile is symmetric and Gaussian 
far from the front, a± = 2. 



V. COMPARISON WITH SIMULATIONS 



To test the two approximations we simulated the PNG 
process. The simulation results presented below are for 
a one-dimensional chain of length L = 10 4 and represent 
a single realization. This study is different than previous 
numerical studies which simulated the PNG process on 
a lattice |24|. Here we treated time and space as con- 
tinuous variables, and this enables comparison with the 
above theories. 

The time dependence of the uncovered fraction for 
the first four layers is shown in Fig. 2. It is seen that 
the MFT and the LRR approaches provide upper and 
lower bounds, respectively, for the actual PNG cover- 
age. Additionally, the MFT provides a better approxi- 
mation for the uncovered fraction, Sj(t). For early times, 
the height and width predicted by either approximations 
are quite close to simulation results, as shown in Fig. 3. 
In fact, both approaches agree to the first significant 
order in time, as both Eq. ( |l2| ) and Eq. ( p3| ) predict 



Sj (t) = 1 — (fjyr^ 2 "' ■ The disagreement between the two is 
of the order i 2j+2 . As the two approximations give upper 
and lower bounds for the PNG process, we conclude that 
this is the leading early time behavior of Sj(t). 




Fig. 2. Uncovered fraction Sj(t) versus t for j — 1,2,3,4. 
MFT and LRR approaches give lower and upper bounds, re- 
spectively, for the uncovered fraction. 




Fig. 3. Short time behavior of the height h(t) and the width 
w(t). MFT is closer to the actual behavior. 

However, both approximations become progressively 
worse at later times. This is due to the fact that the 
asymptotic behavior of the width is predicted incorrectly 
(see Table 1). Fig. 4 shows that at least in one dimen- 
sion, the PNG asymptotic behavior belongs to the KPZ 
universality class. As our simulations are continuous in 
space and time, they enable measurement of the surface 
growth velocity. Although in principle, there is no rea- 
son to expect that the surface growth velocity in an in- 
finite and a finite system are equal, the numerical ve- 
locity Unoncq = L.41 ± 0.01, which corresponds to (non- 
equilibrium) growth on an infinite substrate, is in very 
good agreement with the known analytical equilibrium 
growth velocity v cq = \/2 |2C|] . Numerical values for v, 
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j3 and <r± are summarized in Table 1. Two-dimensional 
simulations give a velocity v% « 1.4 [0] compare to the 
values v 2 = 1.671 (MFT) and v 2 = 1.137 (LRR). This 
suggests that MFT improves for higher dimensions. Fur- 
thermore, the exponent (3 decreases with the dimension 
and it vanishes for d > 4 Above this critical dimen- 
sion, mean-held behavior occurs. 

We also examined the extremal behavior of Sj(t) us- 
ing the simulations. The scaling prediction (||) holds as 
the simulation data is consistent with the exponent value 
er_i_ = 3/2. Given that the PNG uncovered fraction is 
bounded by the two approximations, the parameters v, 
(3, and therefore o+ are similarly bounded (see Table 1.) 
Since the scaling argument involves the width, one can 
not conclude a priory that the same holds for er_ . Nev- 
ertheless, the exponent found in the simulation is quite 
close to 2, or possibly slightly larger. 




10° 10 1 10 2 10 3 



t 

Fig. 4 Long time behavior of the width. Early behavior is 
linear and late behavior is in the KPZ universality class i 1 / 3 . 





MFT 


PNG 


LRR 


Vl 


1.75735 


1.41 ±0.01 


1.12838 


p 





1/3 


1/2 


°+ 


1 


3/2 


2 


<7_ 


2 


> 2 


2 



Table 1: Characteristics of the three approaches for the one- 
dimensional PNG model. 



VI. CONCLUSIONS 

We have investigated a continuum model of multilayer 
growth, the PNG model. We confirmed that the LRR 
approximation implies a dimension-independent asymp- 
totic behavior with the roughness exponent (3 = 1/2. 
Moreover, we found that the full coverage profile ap- 
proaches a dimension independent form, thus implying 
that this approximation ignores interactions. We devel- 
oped an alternative self-consistent mean field approach 
that provides a better approximation for the uncovered 



fraction, and close estimates for the early time behavior. 
Additionally, the two approximate approaches combine 
to give upper and lower bounds for statistical properties 
such as the coverage, the velocity, and the roughness. 

The mean field approach predicts a smooth interface, 
(3 = 0. With increasing the dimension, the roughness 
exponent decreases from /3 = l/2at<i=Oto/3 = Oat 
d> d c , with d c the critical dimension 0]. Thus, for suf- 
ficiently high dimensions the mean-field behavior should 
emerge. 

The above mean-field theory should allow computa- 
tion of space-time correlation functions and structure 
functions. Even more detailed analysis may be possible 
for one-dimensional substrates. A bigger challenge is to 
solve the PNG process analytically. Such a solution will 
undoubtedly illuminate the theoretical understanding of 
non-equilibrium growth. 

PLK acknowledges support from NSF and ARO. 
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